*************************************************************
* Replication codes for analyses in Section 6.2 of
* Beyond Yield Response: Weather Shocks and Crop Abandonment
* May 2020
*************************************************************
	
	use dta\dta_for_irr, clear
		
* Adjustments on variables
		
	* scale units	
		
		replace prec = prec/100
		gen prec2 = prec^2

	* new time var	
		gen t = year-1980
		gen t2 = t^2
		
	* state-year indicator	
		replace stateansi = floor(fips/1000)
		gen stateyear = stateansi*year			
		
	* gen num of yrs for silage
		egen pl_yrs_irr_corn = total(grain_irr_ha_yields!=.), by(fips)
		egen pl_yrs_irr_soy = total(soy_irr_ha_yields!=.), by(fips)
		
	* gen indicators
		gen corn_exist = (grain_irr_ha_yields!=. | grain_irr_ha_acre!=. | grain_irr_pl_acres!=.)
		gen soy_exist = (soy_irr_ha_yields!=. | soy_irr_ha_acre!=. | soy_irr_pl_acres!=.)
		
		replace grain_irr_ha_yields = 0 if grain_irr_ha_yields == . & corn_exist == 1
		replace soy_irr_ha_yields = 0 if soy_irr_ha_yields == . & soy_exist == 1
		replace grain_irr_ha_acres = 0 if grain_irr_ha_acres == . & corn_exist == 1
		replace soy_irr_ha_acres = 0 if soy_irr_ha_acres == . & soy_exist == 1
		
	* taking log
		
		foreach x in grain_irr_ha_acres grain_irr_ha_yields soy_irr_ha_acres soy_irr_ha_yields {
			gen l_`x' = log(`x') 
			}	
		
		gen shr_grain_irr_ha_pl = grain_irr_ha_acres/grain_irr_pl_acres
		gen shr_soy_irr_ha_pl = soy_irr_ha_acres/soy_irr_pl_acres
	
* regressions
		
	estimates clear	
	xtset fips year
	
	* grain spline spec
	
	gen gdd8_26 = (dday8 - dday26)/100
	gen gdd26_35 = dday26 - dday35
		
	qui reghdfe shr_grain_irr_ha_pl c.prec##c.prec gdd8_26 gdd26_35 dday35 ///
		i.stateansi#c.t i.stateansi#c.t2 i.t if pl_yrs_irr_corn>5, ///
		absorb(fips) vce(cluster fips stateyear)		
		estimates store IrCshr_dday_2w	
		
	qui reghdfe l_grain_irr_ha_yields c.prec##c.prec gdd8_26 dday26 ///
		i.stateansi#c.t i.stateansi#c.t2 i.t if pl_yrs_irr_corn>5, ///
		absorb(fips) vce(cluster fips stateyear)		
		estimates store IrCyld_dday_2w	
		
		* irrigated corn estimates (Appendix Table A9)
		esttab IrC*_dday_*, b(4) se(5) keep(*prec* gdd* dday*)			
	
		* effects evaluated at knots for plotting splines (Figure 6)
	
		estimates restore IrCshr_dday_2w
		margins, dydx(prec) at(prec = (0(2)16))
		lincom gdd8_26*(26-8)/100
		lincom gdd8_26*(26-8)/100 + gdd26_35*(35-26)
		lincom gdd8_26*(26-8)/100 + gdd26_35*(35-26) + dday35C*(41-35)
		
		estimates restore IrCyld_dday_2w
		margins, dydx(prec) at(prec = (0(2)16))	
		lincom gdd8_26*(26-8)/100
		lincom gdd8_26*(26-8)/100 + dday26C*(41-26)	
	
	* soy spline sepc

	qui reghdfe shr_soy_irr_ha_pl c.prec##c.prec gdd8_26 gdd26_35 dday35 ///
		i.stateansi#c.t i.stateansi#c.t2 i.t if pl_yrs_irr_soy>5, ///
		absorb(fips) vce(cluster fips stateyear)		
		estimates store IrSshr_dday_2w	
		
	qui reghdfe l_soy_irr_ha_yields c.prec##c.prec gdd8_26 dday26 ///
		i.stateansi#c.t i.stateansi#c.t2 i.t if pl_yrs_irr_soy>5, ///
		absorb(fips) vce(cluster fips stateyear)		
		estimates store IrSyld_dday_2w	

		* irrigated soy estimates (Appendix Table A9)
		esttab IrS*_dday_*, b(4) se(5) keep(*prec* gdd* dday*)		

		* effects evaluated at knots for plotting splines (Figure 6)
		
		estimates restore IrSshr_dday_2w
		margins, dydx(prec) at(prec = (0(2)16))
		lincom gdd8_26*(26-8)/100
		lincom gdd8_26*(26-8)/100 + gdd26_35*(35-26)
		lincom gdd8_26*(26-8)/100 + gdd26_35*(35-26) + dday35C*(41-35)	
		
		estimates restore IrSyld_dday_2w
		margins, dydx(prec) at(prec = (0(2)16))			
		lincom gdd8_26*(26-8)/100
		lincom gdd8_26*(26-8)/100 + dday26C*(41-26)	
		
